Doping dependence of spin and orbital correlations in layered manganites 
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We investigate the interplay between spin and orbital correlations in monolayer and bilayer man- 
ganites using an effective spin-orbital t-J model which treats explicitly the Cg orbital degrees of 
freedom coupled to classical t2g spins. Using finite clusters with periodic boundary conditions, 
the orbital many-body problem is solved by exact diagonalization, either by optimizing spin con- 
figuration at zero temperature, or by using classical Monte-Carlo for the spin subsystem at finite 
temperature. In undoped two-dimensional clusters, a complementary behavior of orbital and spin 
correlations is found — the ferromagnetic spin order coexists with alternating orbital order, while 
the antiferromagnetic spin order, triggered by t2g spin superexchange, coexists with ferro-orbital or- 
der. With finite crystal field term, we introduce a realistic model for Lai_i,Sri+a;Mn04, describing 
a gradual change from predominantly out-of-plane 3z^ — r'^ to in-plane x^ — y^ orbital occupation 
under increasing doping. The present electronic model is sufficient to explain the stability of the 
CE phase in monolayer manganites at doping x — 0.5, and also yields the C-type antiferromagnetic 
phase found in Ndi_a;Sri+i,Mn04 at high doping. Also in bilayer manganites magnetic phases and 
the accompanying orbital order change with increasing doping. Here the model predicts C-AF and 
G-AF phases at high doping x > 0.75, as found experimentally in La2-2iSri+2a:Mn207. 

(Published in Phys. Rev. B 73, 10445i, 2006.) 

PACS numbers: 75.47.Lx, 71.30.+h, 75.10.Lp, 75.75.-|-a 



I. INTRODUCTION 

Colossal magnetoresistance (CMR) manganites are 
characterized by a complex interplay of charge, spin, 
orbital and lattice degrees of freedom. Although man- 
ganese oxides have been known for more than 50 yearsji 
their properties are still not adequately understood. Af- 
ter the observation of the CMR effect f^ the modeling of 
complex behavior in this class of compounds has become 
the focus of intense research activity in the modern con- 
densed matter theoryi^i^ These recent studies demon- 
strate that one has to go beyond the simple double ex- 
change (DE) model of Zener— in order to investigate a 
complex interplay between magnetic, orbital, charge and 
lattice degrees of freedom. 

Similar to perovskite manganites, monolayer 
Lai_2;Sri+3;Mn04 and bilayer La2-2a;Sri+22;Mn207 
manganites have interesting and still poorly understood 
physical properties. The undoped monolayer LaSrMn04 
compound has the same magnetic structure as K2NiF4, 
i.e., it exhibits an G-type antiferromagnetic (G-AF) 
order within the layers i^ This suggests a different or- 
bital state than that realized in three-dimensional (3D) 



LaMnOa perovskite, where A-type antiferromagnetic 
(A-AF) order, with ferromagnetic (FM) ab planes and 
antiferromagnetic (AF) order along c direction, is ob- 
served. Furthermore, unlike in Lai_a;Sra;Mn03, in doped 
monolayer Lai_a;Sri_(_a;Mn04 compounds no FM metallic 
phase was observed, but instead short-range magnetic 
correlations of various types were reported )^i^i^° indicat- 
ing frustrated magnetic interactions. This behavior is 
puzzling and was not explained by the theory so far. 

Also in bilayer La2-2xSri_|_2a;Mn207 manganites a 
competition between magnetic interactions of different 
origin was observed, resulting in rather complex phase 
diagramjiiiiS which is a challenge for the theoretical 
models. At higher doping x ~ 0.45 the magnetic order 
changes from FM to ^-AF phase. The observed phase 
transitions have been ascribed both experimentally^'^ and 
theoreticalljii^ii^ to the varying crystal field splitting be- 
tween Cg orbitals under increasing doping. 

In spite of certain similarities between monolayer and 
bilayer compounds, ^^ the magnetic correlations close to 
half doping {x ^ 0.5) are different. A metallic FM phase 
is observed in bilayer manganites up to a; '--^ 0.45,^^ while 
it is absent in monolayer compounds. While the CE-type 



AF order is quite pronounced at a; ~ 0.5 in monolayer 
manganitesji^ the A-AF phase is instead more stable in 
bilayer systemsfii*i^ and the charge order and orbital bi- 
stripes were also observed for higher doping 0.55 < a:: < 
0.6.19 

In the present paper, we intend to focus on the role of 
orbital degrees of freedom in stabilizing various types of 
magnetic order observed in monolayer and bilayer man- 
ganites. The behavior of Cg electrons is dominated by 
large Coulomb interaction U. Therefore, we employ an 
effective spin-orbital t-J model similar to that derived 
for the one-dimensional (ID) chain)2£ and generalize it to 
the present situations. Thereby, we implement also the 
Hund's exchange interaction Jh which enforces the spin 
s = 1/2 of an Cg electron to follow the t2g spin S = 3/2 
at each site in the ground state. Unlike in the ID case, 
the orbital Cg flavor is not conservedfSl which enhances 
quantum fluctuations. They contribute to intersite cor- 
relations and we show that a close relationship between 
orbital and spin correlations nevertheless persists in the 
ground state and at low temperature. 

The previous theoretical studies revealed a competi- 
tion between different types of magnetic order. Among 
them the most spectacular ones are phases with coex- 
isting FM and AF bonds: (i) E-type AF phase in un- 
doped systemSfSSiS^ and (ii) the CE phase at half doping 
{x = 0.5). The former one has been experimentally ob- 
served for the very strongly Jahn- Teller (JT) distorted 
case only, which is at variance with its theoretical pre- 
diction for undistorted compounds. The mechanism of 
stability of the CE phase is also still under debate. While 
it has been shown that local JT distortions induce the 
CE-type AF oidei^^i^ it remains unclear whether it 
could follow from electronic interactions alone. It was 
argued before that this complex type of magnetic and 
orbital order (00) might either originate from conflict- 
ing phasesr& or could be stabilized by intersite Coulomb 
interactions^ Here we address these various mechanisms 
proposed before by presenting the evidence obtained by 
numerical simulations of flnite clusters within a realistic 
electronic model including Coulomb interaction. We also 
show that this model gives a satisfactory description of 
magnetic correlations over the entire doping range. 

The paper is organized as follows: In Sec. II we present 
the effective t-J orbital model in the regime of large U 
for Cg electrons, moving in two-dimensional (2D) clus- 
ters simulating monolayer manganites, or in \/8 x \/8 x 2 
clusters standing for bilayer manganites. We also present 
shortly two numerical methods: the exact diagonaliza- 
tion with Lanczos algorithm used to solve the orbital 
model for fixed spin conflgurations at zero temperature 
(T = 0), and its combination with Monte Carlo simu- 
lations of spin core (i2g) conflgurations, which leads to 
a coupled spin-orbital problem at T > 0. In Sec. Ill 
the model for monolayer manganites is analyzed in dif- 
ferent doping regimes. We report the phase diagrams 
obtained for undoped and half doped systems, and re- 
late the obtained magnetic phases to orbital occupations 



and intersite orbital correlations. Thereby we highlight 
the interrelation between spin and orbital order and their 
dependence on increasing doping. The study of bilayer 
manganites in Sec. IV is limited by the size of Hilbert 
space for the smallest relevant \/8 x \/8 x 2 clusters, so 
we discuss only undoped (x = 0), half doped {x = 0.5) 
and strongly doped (x > 0.8) systems. Finally, in Sec. 
V we summarize the numerical results and present gen- 
eral conclusions deduced from the present study for the 
physical mechanisms operating in layered manganites. 



II. THE MODEL AND NUMERICAL METHODS 

A. Orbital t-J Model 

The effective orbital t-J model described below follows 
from the model of interacting eg electrons. 
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The form of the kinetic energy H^. depends on the se- 
lected basis of orthogonal orbitals, as discussed in detail 
by Feiner and 01es.-2i Here we use the conventional basis 
which consists of 
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orbitals. In a 3D (or bilayer) manganite the kinetic en- 
ergy takes the form. 
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and the last term is absent for a monolayer. The largest 
hopping element t stands for an effective (dda) element 
between two \z) orbitals along the c axis, and originates 
from two consecutive d-p transitions over oxygen orbital 
between the neighboring Mn ions. The same hopping el- 
ement t couples two equivalent directional orbitals along 
either a or 6 axis, i.e., 3x^ — r^ or 3j/^ — r^ orbitalSiSi 

By considering the structural data, a uniform crys- 
tal field splitting of Cg orbitals is expected both for a 
monolayer (^ and for a bilayer systemii^ Therefore, we 
introduce the term. 
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where niaa — claa'^iaa ^^ the electron number operator in 
a = X, z orbital with spin a at site i. This term describes 
the crystal field splitting of eg orbitals which follows from 
the geometry of layered manganites and removes the or- 
bital degeneracy. If i^z > 0, as in undoped LaSrMn04f^ 
the \z) orbitals are favored. 



The electron interactions between Cg electrons are 



given by, 
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The spin operators, Sia = {s^^, Sj^, sf^}, are defined by 
fermion operators in a standard way. 
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The on-site interactions in iJjnt are rotationally invariant 
in the orbital space and are thus described with only 
two parameters!^ the Coulomb element U, and a Hund's 
exchange element Jh- The intcrsite interactions oc V, 
where each bond (ij) is included only once, do not depend 
on the orbital type and thus involve total electron density 
operators, nt = X^aij^^ao-- 

The Hamiltonian ^ does not include yet the t2g elec- 
trons which could in principle be described again by 
a similar multiband Hamiltonian.'^^ However, in reality 
t2g electrons localize due to large Coulomb interaction 
U ^ t, so it suffices to consider local spins S = 3/2 
built by three t2g electron spins of either Mn'^+ or Mn'*+ 
ion. In both cases each t2g orbital is singly occupied and 
virtual intersite hopping processes contribute to the su- 
perexchange oc J'. In this parameter regime, we also 
find hopping processes of Cg electrons between two Mn'^"'" 
ions, dfd'^ ^ d^d^, that lead to e^ configurations at site 
i (with e'i electrons cither in one or at two different or- 
bitals) and generate the Cg part of the superexchange 
oc J = it'^/U. When the system is doped, direct hop- 
ping processes are also possible for Mn'^+-Mn''^+ pairs, so 
one finds an effective t-J model, similar in spirit to the 
spin t-J model derived almost three decades ago from the 
Hubbard modelj^S 

An important difference to the Hubbard model, how- 
ever, arises due to Hund's exchange, —iJnSi ■ Si which 
favors high spin states at each site. Hund's exchange 
is closed to the atomic value of Jr ^ 0.9 eV, so is 
sufficiently larger than the hopping, estimated for the 
bilayer— to be t ~ 0.48 eV, to assume that spins of Cg 
electrons are aligned with the S = 3/2 core spin formed 
by the t2g electrons. Therefore, core spin determines a 
local spin quantization axis at each site i, the spins of 
Cg electrons can be integrated outj2Z. and the effective 
orbital t-J model takes the form. 
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This Hamiltonian may be also obtained by generalizing 
the effective ID model of Ref. |23 to doped layered man- 
ganites. It depends not only on the microscopic param- 



eters introduced above, but also on the actual configu- 
ration S of t2g spins on the lattice which determine the 
hopping term Ht by the double exchange mechanism. 

The first term Ht in Eq. {TJ stands for the hopping of 
Cg electrons in the limit of large on-site Coulomb inter- 
action U ^ t. Its form is given below for the monolayer 
and bilayer system separately. We emphasize that by 
performing rigorous projection to the C/ — > cx) limit'^^ the 
orbital degree of freedom of Cg electrons survives, but in 
agreement with the basic idea of the spin t-J model^S 
and in contrast to Eq. Q, the hopping processes are 
now limited to the subspace without (intraorbital and in- 
terorbital) double occupancies. The double occupancies 
generated in virtual charge excitations by either eg or t2g 
electrons contribute in second order of the perturbation 
theory and give the superexchange interactions ex J or 
oc J', respectively, contained in Hj and Hj' terms. 

Similar to LaMnOs,"^^ the superexchange in undoped 
LaSrMn04 and in La2-2a;Sri+2a;Mn207 is given by a su- 
perposition of several terms. The Cg term Hj favors ei- 
ther FM or AF spin order on a bond (ij), depending on 
the pair of occupied Cg orbitals at neighboring sites i and 
j. This makes the form of the superexchange term Hj 
depend both on the actual geometry in layered mangan- 
ites, and on the used orbital basis. As the hopping term 
Ht, its form depends on the system and is given below. 

In contrast, the superexchange interactions induced 
by charge excitations of t2g electrons Hj' are identi- 
cal for planar and bilayer manganites. They are fre- 
quently treated as an effective AF superexchange be- 
tween S — 3/2 core spins, although they couple de facto 
two manganese ions in high-spin configurations, and thus 
depend on the actual total number of d electrons at both 
ions.'^'*- We have verified, however, that the t2g superex- 
change terms derived for these different configurations 
are of the same order of magnitude, so in a good ap- 
proximation one may indeed simulate their effect by the 
Heisenberg Hamiltonian with an average exchange con- 
stant J' > between S = 3/2 core spins which favors 
AF spin order; this interaction is described by the term. 
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For convenience, we use classical core spins Si of unit 
length (compensating their physical value S = 3/2 by a 
proper increase of J'), i.e., we replace the scalar products 
of spin operators on each bond by their average values. 
The classically treated Si are represented by polar angles 
{'di, (f)i} — then the spin product on a bond (ij) is given 
by: 



{S,-S,)^S^{2\u.,,\^-l) 
where the spin orientation enters via 
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depending on the angle Oij between the two involved 
spins and on the complex phase Xij ■ 

The remaining terms in Eq. 0, H^ and Hy, stand 
for the crystal field splitting of two Cg orbitals caused by 
geometry and for the nearest neighbor Coulomb interac- 
tion. The first term, 
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involves projected density operators, fii^ and n^j,, that 
act in the restricted Hilbert space without (interorbital 
and intraorbital) double occupancies. The latter term, 
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(12) 



is the intersite Coulomb repulsion which survives as the 
only term from Eq. ^. ft plays a role at finite doping 
where it can induce charge order. Here fii is the total 
Cg electron density operator at site i which acts in the 
restricted Hilbert space: 

fii = n„ +ni:,. (13) 

By construction (hi) < 1. We include this term in half 
doped manganites and investigate in Sees. If f f Cl and lHI Dl 
whether it helps to stabilize the CE phase with coexisting 
charge order. 



B. Monolayer manganites 

We start with the general form of the hopping term 
Ht which follows from Eq. Q for planar (monolayer) 
manganites in the limit of [/ ^ i. 
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where an operator c^^ creates an electron in \a) state at 
site i when it is unoccupied by any other (a or j3) electron, 
i.e., implements rigorously the restriction of the Hilbert 
space to the subspace with no double occupancies. Fur- 
thermore, Hund's exchange between the itinerant Cg elec- 
trons to the t2g core spins is large enough to justify the 
restriction to the subspace of tg electrons parallel to the 
local t2g spin, i.e., we treat spinless fermions with an or- 
bital flavor a = x, z, see Eq. ^. Electrons with antipar- 
allel spins and double occupancies are treated in second 
order perturbation theory, see below. 

In agreement with double exchange mechanismjSiS^ the 
effective hopping strength is modulated by the scalar 
product of the core spins at the respective sites — it 
is maximal for parallel spins and vanishes when spins are 
antiparallel. Accepting the classical treatment of inter- 
site correlations between core spins, we use 



with Uij given by Eq. H10|l . The first factor tap is the 
orbital-dependent hopping strength. Its form depends 
on the used orbital basis, and for a given pair of orbitals 
{a,/3} depends on the direction of the bond {ij), as ex- 
plained below. 

While a representation using three directional orbitals 
IC) along each cubic direction gives the simplest expres- 
sion for the hopping of eg electrons, ^"'^ it is more conve- 
nient to consider here the fixed orthogonal orbital basis 
given by Eq. (0) in a 2D geometry, for which the orbital- 
dependent hopping strength tap in Eq. H15|) is given by: 
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for a bond along the a and h cubic axis, respectively. The 
hopping term takes then the form. 
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The correlated fermion operators c]^ — '^Izi^ ~ ''^ix) f^nd 
c'-^ = cj^(l — riiz) act in the restricted Hilbert space 
and create a \z) {\x)) electron only when site i is ini- 
tially empty. The hopping term H17|l describes therefore 
spinless fermions with an orbital flavor in a restricted 
Hilbert space, in analogy to the original t-J model in spin 
space,''*' but with an anisotropic hopping term. We em- 
phasize that the orbital flavor is not conserved along the 
hopping process in dimensions higher than onei^ This 
has important consequences for the phase diagram of the 
orbital Hubbard model, and destabilizes the 00 in doped 
manganites. ^^ 

Similar as for the undoped LaMnOs (see R ef . l3a) . the 
Cg superexchange term for undoped LaSrMn04 depends 
on the pair of occupied Cg orbitals at sites i and j for 
a given bond (ij). If the spins are in the FM configu- 
ration, the superexchange interactions do not vanish but 
are simply reduced to purely orbital interactions which 
favor alternating directional (3^^ — r^-like) and planar 
(x^ — y^-like) orbitals along every cubic directioni2I An 
effective orbital superexchange model presented below 
originates from the complete spin-orbital model^ that 
included the complete multiplet structure of the excited 
d*5 and d * 4 states;^ and focuses on the orbital dynam- 
ics in the presence of spin fluctuations which influence 
orbital superexchange interactions. 

The superexchange due to Cg electron excitations con- 
tains spin scalar products multiplied by orbital inter- 
actions on the bondS)^ and the full many-body prob- 
lem would require treating the coupled spin and orbital 
dynamics. Here we study only the orbital correlations 
and their consequences for the magnetic order by replac- 
ing the scalar products of spin operators on each bond 
by their average values ©, as done before in the ID 



model, ^^ which is equivalent to decoupling spins and or- 
bitals in a mean- field approximation. As a result, one is 
treating the orbital many-body problem coupled to the 
classical spins. 

The spin operators are replaced by their expectation 
values following Eq. H1Q|I . In the present case of a mono- 
layer one finds the Cg superexchange term. 
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where number operator tuq refers in each case to the di- 
rectional orbital \C} along a given bond (ij), and 
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depend on the bond direction, with the sign — (-I-) in 
Eq. H19|l corresponding to a (b) axis. The operators are 
defined by orbital T — 1/2 pseudospin operators. 
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with two eigenstates of Tf , see Eq. (01. The superex- 
change constant J = t^ /e{^Ai) is determined by the low- 
est high-spin excitation energy e(M.i) at the Mn^+ ioni2^ 



C. Bilayer manganites 

The hopping term for the bilayer manganites can be 
written again in terms of directional orbital along the 
bonds, by extending Eq. H17|) . In fact, the interlayer 
hopping term along c-axis is then diagonal in {jo:), |z)} 
basis. 
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As for a single plane, the hopping Ht describes spinless 
fermions with an orbital flavor in a restricted Hilbert 
space. In each plane the orbital flavor is not conserved 
along the hopping processiSi 

We apply again the same approach described above to 
the superexchange interaction; it leads in the present case 



of a bilayer system to the following expression for the eg 
superexchange. 
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Compared to Sec. IIIBI it is extended by the interlayer 
coupling terms, with the charge and orbital operators 
given by fiiz and Tf along c axis [the remaining terms 
for the bonds in ab planes are the same as in Eq. (|18l) ]. 
For simplicity, we neglect here the difference between the 
intralayer and interlayer hopping elements, although for 
a detailed comparison with experiment it might be nec- 
essary to include different bond lengths for each cubic 
direction. 



D. Density distribution and intersite correlations 

Below we present the physical quantities which are 
used to characterize electron density distribution and in- 
tersite spin and orbital correlations in the ground states 
and at low temperature. The density distribution in both 
monolayer and bilayer systems at doping x is described 
by average electron densities in orbital a, 
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with N being the number of lattice sites. 

Depending on the actual parameters, we encountered 
several different types of magnetic order, both in mono- 
layer and in bilayer model. The FM and AF spin interac- 
tions compete with each other, so one of them is selected 
by a particular type of orbital correlations. The possi- 
ble generic types of magnetic order, all observed in the 
manganites j2i^ are shown schematically in Fig. ^ We 
have investigated intersite correlations at short distance 
r by: (i) intersite spin correlations, given by the scalar 
product lO of the classical core spins. 
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and (ii) intersite orbital correlations. 
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FIG. 1: (Color online) Schematic spin structures in the mag- 
netic phases of doped manganites. The number of FM bonds 
decreases gradually from FM to G-AF phase, through A-AF 
and C-AF phase, with eight and four FM bonds in a cube. 
For an ab monolayer, FM and ^-AF phases are equivalent. 



The orbital operators are defined using a particular or- 
bital basis, 



T,{0) ^T^ cos 61 -I- Tf sine. 
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and the pscudospin operators {T^ ,T'f} at site i are given 
by Eqs. H20() and (|21ll . For instance, 9 — corresponds 
to T^^Tf^^, = 7r/2 - to T-T-^^, and 9 = 2^3 - to 
T/T^, -, with C standing for the directional orbital along a 
axis. The orbital correlations expected in undoped man- 
ganites are of AO type on two sublattices, which suggests 
that the orbital correlations Tg{'r) defined as in Eq. H27|l 
are predominantly negative for nearest neighbors. These 
correlations were investigated along the (10) and (11) di- 
rections in 2D clusters. In the bilayers we will distinguish 
between two different cases for both above directions, and 
consider: either (i) both sites within the same ab layer, 
or (ii) two sites which belong to different layers. 

Spin correlations in a doped system were uniquely de- 
termined by core spin correlations. When Cg holes are 
present, it is also of interest to investigate spin correla- 
tions near the hole, so we show in some cases the corre- 
lations between the two nearest neighbor spins separated 
by a hole, 
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for bilayer ones. Here, the spin configurations were fixed, 
corresponding to different possible magnetic phases: FM, 
AF, C-AF, £'-AF, and CE phase. To solve the orbital 
problem specified by the selected core spin configuration, 
we employed the Lanczos algorithm, which is suitable for 
very large matrices treated here and guarantees fast con- 
vergence to the ground state in each case. We then deter- 
mined the global ground state by comparing the ground- 
state energies obtained for different magnetic phases. 

At finite temperatures, we investigated the effective 
orbital t-J model by making use of a combination of 
Markov chain Monte Carlo algorithm for the core spinsSS 
with Lanczos diagonalization. For each classical core spin 
configuration occurring in the MC runs, we defined the 
actual values of classical variables {wjj}, and next em- 
ployed Lanczos diagonalization to solve the many-body 
problem posed by the orbital model. In each case we 
obtained the free energy for that core spin configura- 
tion from the few lowest eigenstates, which was next 
used to decide acceptance in the MC runs. We mea- 
sured the Boltzmann weight of these lowest eigenstates 
and thereby made sure that only states with negligible 
weight were discarded. We emphasize that a complete 
many-body problem in the orbital subspace was solved, 
so our method differs from the standard algorithm used 
for noninteracting electrons which employs free-fermion 
formulae, as e.g. in Ref. i3^ 

In the MC updates, the angle of spin rotation was opti- 
mized to keep acceptance high enough. If acceptance was 
very good, several rotations were performed in each up- 
date. From time to time, a complete spin flip S^ -^ — S; 
was proposed. Autocorrelation analysis was employed 
to obtain reliable error estimates and several hundred 
effectively uncorrelated samples were considered, taking 
particular care of burn-in and thermalization processes. 
Finally, wherever autocorrelations were observed to be 
particularly long, e.g. in symmetry-broken states like 
the CE phase, the method of parallel tempering'^^ was 
employed. 



where (1 — fn) stands for the hole density at the central 
site i, and two spins Si±s occupy two adjacent sites in 
either a or b direction. 



III. NUMERICAL RESULTS FOR 
MONOLAYER MANGANITES 

A. Undoped 2D clusters 



E. Numerical methods 

We employed two different but related numerical meth- 
ods to investigate finite clusters described by the orbital 
t-J model: (i) exact diagonalization at zero temperature 
(T = 0), and (ii) its combination with Markov chain 
Monte Carlo (MC) at finite temperature T > 0. The 
ground state at T = for representative electron fill- 
ings was determined by solving the orbital t-J model for 
several possible types of spin order, using clusters with 
periodic boundary conditions: \/8 x -\/8 and 4x4 clus- 
ters for monolayer manganites, and \/8 x \/8 x 2 clusters 



In undoped manganites (at a; = 0), the hopping is en- 
tirely blocked by large Coulomb interaction U, and the 
nearest-neighbor Coulomb interaction V H12(l is irrelevant 
as the electron density is n = 1 at each site. Therefore, 
the Hamiltonian for a monolayer Eq. Q reduces to the 
superexchange terms given by Eqs. (jHl and p8|l . accom- 
panied by the crystal-field splitting of eg orbitals, given 
byEq. ((TT|l . Therefore, the ground state for finite J is de- 
termined by only two parameters: J' /J and Ez/J- The 
phase diagram of the spin-orbital model (0 in {J',Ez) 
plane obtained by exact diagonalization of a 4 x 4 cluster 
at r = is shown in Fig. |21 
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FIG. 2: Phase diagram for the undoped layered manganites 
(a; = 0), with the stabihty regions of FM, G-AF, and C-AF 
phases (see Fig. in {J',Ez) plane, as obtained at T = 
with a 4 X 4 cluster. The core spin superexchange J' and the 
crystal field splitting Ez axe given both in units of J, and in 
units oft for J = 0.125t. 




FIG. 3; (Color online) Orbital correlations Te{r) ^ in 
(01) direction in the FM phase, as obtained for the undoped 
4x4 cluster using T/ operators defined for different bases 
of orthogonal orbitals: ^ = — {|a;),|z)}; 6 = 7r/2 — 
{{\x) ± \z))/V2}; e = 27r/3 — {10,10}- The latter choice 
involves directional orbitals within the ah plane. Orbital cor- 
relations shown in (11) direction for 6 — ■k/2 demonstrate an 
almost perfect AO order. Parameters: Ez = J' = 0. Spin 
correlations obtained for the spin model (i.e. the Heisenberg 
model for s = 1/2) on a 4 x 4 cluster at n = 1 are shown for 
comparison by squares and diamonds. 



A particularly simple result is obtained at J' = -B^ = 0, 
where the FM phase has the lowest energy. As Uij = 1 
in a ferromagnet, only the first term in the eg superex- 
change H18|l survives and stabihzes the alternating orbital 
(AO) order. Of course, this phase is also stable in an un- 
physical regime^ of negative J' < in a broad range of 
crystal field splitting, see Fig. [3 This coexistence of the 
AO order with the FM spin correlations is generic — it 
confirms the trend observed before in the ID model^^ and 
on ladderSjii and agrees with the Goodenough-Kanamori 
rulesi^ When robust AO state develops at the orbital de- 
generacy {Ez = 0), the FM phase extends to a broader 
range of J' > than at \Ez\ > 0. 

As the FM and AF terms in eg superexchange H18|) 
compete with each other, the G-AF phase has only 
slightly higher energy than the FM one at Ez = 0, so 
it can be rather easily stabilized by finite AF t2g su- 
perexchange J' > 0. One finds indeed a transition to 
the AF order at J' > 0.009U (we use below t = 1 as 
an overall energy unit). Note, however, that near orbital 
degeneracy {E^ ~ 0) robust AO state develops as both 
orbitals can participate in equal amount, and therefore 
the FM phase extends here to higher values of J' than 
at large orbital splitting \Ez\. Therefore, the crystal field 
energy may easily tip the energy balance and stabilize 
a different type of magnetic order, although Ez controls 
primarily the charge distribution. Indeed, at large \Ez\ 
only one of eg orbitals is selected, AO order is hindered, 
so the AF phase is favored and occurs already for smaller 
J' than at Ez = 0, see Fig. |21 For Ez < 0, \x) orbitals 
are occupied and one finds large AF eg superexchange 
— therefore in this case the AF phase is favored even 



for J' = 0. In contrast, when \z) orbitals are favored 
for Ez > 0, the eg superexchange is weaker by a factor 
of {txx/tzzY — 9, and the G-AF phase can be stabilized 
only by J' ~ O.Oli. Moreover, one finds that the com- 
petition between the FM and AF order may induce the 
G-AF phase in the crossover regime between the FM and 
AF phase. We note that the G-AF phase is expected to 
be further stabilized by finite cooperative JT distortions. 

In order to get more insight into the phases of Fig.|21 we 
investigated orbital correlations between first and second 
neighbors. Let us first look at the ferromagnetic phase at 
Ez = 0: Since there is no hopping at a; = and since all 
bond are FM (i.e. Uij = 1), the Hamiltonian comprises 
only the first term of Eq. ^1 While it looks similar to 
a spin t-J model (or in this undoped case the Heisen- 
berg model for s = 1/2), the important difference is that 
the orbital model is not isotropic. As a consequence, the 
orbital correlations H27|) differ markedly from those fa- 
miliar from the isotropic spin t-J model, and depend on 
the considered orbitals basis. They are shown in Fig. O 
for various orbitals parameterized by angle 9 [sec Eq. 
(|28|l ]. and one finds that almost no correlations are found 
when angle 61 = is selected, i.e., for 7^j(0) = (T^^T/). 
In contrast, orbital correlations designed to detect the 
AO order, with angle 9 ~ 7r/2, are quite distinct, as for 
instance for 9 = 27r/3, when correlations between two 
directional ix^ — r^ /iy^ — r^ orbitals within the plane 
are measured by 7ij(27r/3) — T-'TJ. As this correlation 
function is negative, one finds the 00 close to staggered 
directional orbitals on two sublattices, with 9i^A = 9 and 
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FIG. 4: (Color online) Spin correlations S{f) 1261 1 in the un- 
doped monolayer, as obtained in MC simulations with a 4 x 4 
cluster for J' — and E^ — (FM correlations), and for 
J' =0, E^ = -0.5i, as weU as for J' = 0.05i, E^ = (both 
sets give AF correlations). Error bars are smaller than symbol 
sizes. Parameter: J = 0.125f, f3t = 100. 
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Remarkably, the orbital correlations defined by the di- 
rectional orbitals on two sublattices (—0.1816) are very 
similar to those obtained with 7ij(27r/3) (-0.1892), but 
both types of 00 do not minimize the negative orbital 
superexchange. A common large negative contribution 
to the above values comes from {T[T^) = —0.1854 cor- 
relation — it decreases further when the angle 9 is varied 
towards 9 = 7r/2, where the orbital correlations reach a 
minimal value, %j{'k/2) = {TfT^) = -0.2472. Thus, 
in agreement with earlier findings r^di the quantum cor- 
rection to the classical value —0.25 is very small indeed 
due to the gap which opens in orbital excitations in the 
present 2D case. This means that robust AO order with 
orbitals of the form (|a;) ± \z))/y/2 is realized in an un- 
doped monolayer without crystal field splitting, similar 
to the 00 in the 3D modelfS if this monolayer has FM 
spin order. 

For comparison, we also included in Fig.|31the spin cor- 
relation H2()|) of the spin i- J-model — it is isotropic due 
to SU(2) symmetry, and quantum fluctuations keep the 
value of 5(1) for (10) direction well above the classical 
value —1/4. In contrast, for the orbital model the optimal 
correlations found for 9 = 7r/2 almost attain the classical 
value and thus show almost perfect 00. An additional 
advantage of this robust 00 for the present calculations 
is that the ground-state energy Eq hardly depends on 



FIG. 5: (Color online) Spin correlations S{r) Ij26|l in the un- 
doped monolayer as in Fig. 21 but for decreasing J' from 
J' = (AF correlations) to J' = -0.04f (FM correlations). 
Parameters: J = 0.125t, E^, = -0.5t, pt = 100. 



the cluster size 
at E:, = 0: ^0 = 



one finds for the FM phase {uij = 1) 
-0.21970i for ^/8 x VS, -0.21968t for 



^lOx VlO, and -0.21967t for 4x4 cluster, i.e., finite-size 
effects are negligible. When the AF order is considered 
instead {uij = 0), the second term in the orbital superex- 
change Eq. H18|l dominates and thus almost only in-plane 
\x) orbitals are occupied. Finite size effects are here again 
smafl, with: Eq = -0.18335t for a/8 x ^8, -0.18334i for 
\/lO X VlO, and -0.18333i for 4 x 4 cluster. 

In the present model one does not find the £^-AF phase, 
which has been experimentally observed for the strongly 
JT distorted HoMnOs."*^ We attribute this to the fact 
that the present model does not include phonons and 
thus rather describes manganites with little or no JT 
distortion. The E-KY phase was found before in MC 
simulations using a model neglecting on-site Coulomb 
repulsion in the regime of large J' and small electron- 
phonon coupling A, where it was stabilized by the kinetic 
energy,^'^'^'* It could argued, however, that this is not a 
realistic description, as local Coulomb repulsion (intraor- 
bital U and interorbital U') inhibit Cg electron motion for 
X — Q {n = 1) limit, so the microscopic mechanism of E- 
AF phase in HoMnOa remains puzzling. 

While ground state calculations comparing different 
ordered phases lead to valuable insights at relatively low 
computational cost, one has to realize that such cal- 
culations at T = are still necessarily limited by the 
authors' imagination concerning the magnetic phases to 
consider. We therefore complemented them by unbiased 
MC simulations for a few parameter sets at low temper- 
ature (3t — 100. Figure 01 shows the spin correlations 
(for core t2g spins) obtained from the MC data for one 
[E^ = J' = 0) example of the FM, and two {E^_ = 0, 
J' = 0.05t and E,^ = -0.5i, J' = 0) of the AF phase. 
These data are complemented by the orbital correlations 
(not shown) which correspond closely to the ground state 
results. Both spin and orbital correlations weaken with 
rising temperature, as discussed elsewhere,^^ but for the 
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FIG. 6: (Color online) Electron density in the \z) orbital per- 
pendicular to an ab plane of a monolayer manganite for in- 
creasing crystal field splitting E^, see Eq. IIIH . as obtained 
from MC simulations of undoped (n = 1) \/8 x \/8 clusters 
with various values of J'. Parameters: J = 0.125t, /3t = 100. 



realistic values of J' ~ 0.02t spin correlations melt some- 
what faster even in the present case when the orbital 
interactions induced by the JT effect are neglected. 

The situation is somewhat different for the C-AF phase 
obtained by exact diagonalization at T = for E^ < 
—0.2t (see Fig. |2l. We analyzed the spin correlations Sij 
at low temperature for Ez = — 0.5t and a few selected 
values of J' = 0, -0.022i, -0.03t, -OMt (Fig. \^. For 
J' = one finds AF correlations and for J' = -0.04i FM 
ones, both in accordance with the ground state phase di- 
agram of Fig. 121 However, the phase transition between 
these phases (the value J' — — 0.022i lies right in the 
middle of the C-AF region) does not occur via the C- 
AF phase at finite temperature, which would give an AF 
signal for r = (1, !)■ Neither do the intermediate values 
of J' exhibit the _E-AF phase, which should show an AF 
signal for r — (0, 2) and has only slightly higher energy 
than the C-AF phase at T = 0, nor do the MC snap- 
shots for this parameter range and /3t — 100 suggest any 
other ordered phase. However, the temperature (3t — 100 
might be still too high to allow for longer range magnetic 
correlations precisely in the crossover regime. In conclu- 
sion, we have established that the AF and FM phases of 
Fig. 12 are well supported by the MC results, while we 
believe the transition between them might occur either 
via the C-AF phase, or with phase separationc2i 

In order to understand better the phase diagram of 
Fig. 121 it is instructive to consider the orbital occupation. 
When positive crystal field {E^ > 0) is applied, \z) or- 
bitals which stick out of the ab plane are favored, see 
Fig. El For the AF phase at J' = O.Obt, relatively small 
value Ez = 0.15< suffices already to switch the density 
distribution from mainly |a;) (at E^ — 0) to mainly \z) 
orbital occupation. The width of the crossover regime 
from I a;) to \z) orbital occupation increases with decreas- 
ing J' . As discussed above, the FM phase is characterized 



by an approximately equal occupation of both orbitals 
at Ez = 0. When both Cg orbitals mix, the AO order 
found in the FM phase at J' = is robust and a rather 
large value of crystal field Ez ~ OAt is needed to enforce 
complete \z) polarization. This large value of Ez refiects 
the cooperative character of the FM spin and AO order. 
At the same time, a negative crystal field Ez ^ — O.li 
is needed in order to enforce complete \x) polarization, 
while this transition occurs a.t Ez — for J' — O.Q2t. 

Note that in the ferro orbital (FO) phase with \x) or- 
bitals occupied (at Ez < ~0.2t, the AF spin order is 
induced even for J' = 0, see the spin correlations in 
Fig. ^ and the phase diagram of Fig. (21 This provides 
another example of complementary spin and orbital cor- 
relations that agree with Goodenough-Kanamori rules. 4S 
At J' = 0.02i one finds an intermediate situation in this 
respect — while both at Ez < and Ez > 0.25t the 
the magnetic order is AF and dictated by J', with ei- 
ther \x) or \z) orbitals being filled, in the transition re- 
gion a.tO<Ez< 0.25i one finds AO {\x) ± \z))/V2 or- 
der, but rather unclear magnetic structure at /3t — 100. 
This shows again that magnetic correlations are typi- 
cally lost at lower temperature, particularly when differ- 
ent conflicting trends in magnetic interactions compete 
with each other. 

Experimentally, the AF order with eg electrons oc- 
cupying mainly \z) orbitals was found^ in a monolayer 
undoped compound LaSrMn04, which we interpret as 
a consequence of sufficiently large positive crystal field, 
Ez ^ 0.5i, induced by the 2D structure. With growing 
temperature \x) occupation riseS)^ as we also observed in 
the MC simulationsi^ 



B. Stability of the CE phase at half doping 

Next we investigate the spin correlations at half dop- 
ing (x — 0.5) and determine the range of stabil- 
ity of the CE phase, which was observed in a mono- 
layer Lao.5Sri.5Mn04 compoundAi^ We have performed 
ground state (T — 0) calculations using three clusters 
shown in Fig. [71 with periodic boundary conditions. Be- 
cause of the large Hilbert space at half-filling, MC sim- 
ulations could be completed only for -s/S x -/S clusters 
[for the cluster geometry and one possible realization of 
the CE phase see Fig. [7fa)]. In these simulations, we 
obtained the CE phase for some parameters, e.g. for 
J = 0.125i, Ez = 0, either for J' = 0.05i and V = 0, 
or for J' = 0.025i and V = t. A typical MC snapshot 
is shown in Fig. |S1 As we have used periodic boundary 
conditions, one recognized the CE phase in an \/8 x VS 
cluster repeated several times. 

The spin correlations resulting from the MC runs are 
compared to those of the ideal CE phase (exact diago- 
nalization at T = 0) in Fig. j^l and one finds an almost 
perfect agreement — in both cases, AF and FM signals 
cancel along (01) direction, because of the different pos- 
sible realizations of the CE phase, but along the diagonal 
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FIG. 7: (Color online) Clusters used to investigate the CE 
phase at half doping (x — 0.5): (a) a/8 x a/8, (b) 4x4, and 
(c) 2x8. Dark (light) shading indicates possible realizations 
of staggered FM zig-zag chains in the CE phase. 



(11) direction the spin correlation reaches —0.5. (For the 
Monte Carlo simulations in these parameter regimes, we 
used parallel tempering'^^ in order to sample the different 
symmetry-related realizations correctly.) 

On the one hand, the occurrence of the CE phase 
for classical spins even without the cooperative JT ef- 
fect in the present study is in contrast to the results ob- 
tained recently for quantum mechanical S = 1/2 core 
spins^ where quantum fluctuations suppress the CE 
phase. While it is not completely clear whether S = 3/2 
core spins are better approximated by more quantum 
S — 1/2 or by classical S* ^ oo spins, the latter has 
been shown to be an excellent approximation in the one- 




FIG. 8: (Color online) MC snapshot obtained for doping x — 
0.5, with spin directions indicated by lines. Four y/S x a/8 
clusters are shown (heavy lines indicate the original cluster) 
and one clearly recognizes the magnetic order in the CE phase. 
Parameters: J = 0.125i, J' = 0.05i, E^, = V = 0, I3t = 100. 
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FIG. 9: (Color online) Spin-spin correlation S{r) i'liil 12511 
for the ideal CE phase (CE) and the MC data (MC) for an 
a/Sx a/8 cluster. There is one AF signal at (1, 1). Parameters: 
J = 0.125t, J' = 0.05t, E^ = V = 0, pt = 100. 



orbital model in one dimension. ^^ On the other hand, a 
similar model, but including phonons and without on-site 
Coulomb repulsion between the Eg electrons was treated 
on a 4 X 4 site cluster [see Fig. [7Jb)] in Ref. [5^ Be- 
cause the on-site Coulomb repulsion was not included, 
unrealistically large J' > 0.2t or rather strong electron- 
phonon coupling A > 1.75 were necessary to stabilize the 
CE phase. Although this result could still be improved 
by considering larger clusters, it highlights already the 
importance of strong Coulomb interactions for obtaining 
physically relevant results. 

When FM and AF interactions occur simultaneously 
at different bonds in ab planes, one has to investigate 
the stability of the CE phase by comparing it with the 
C-AF phase which has the same amount of FM and AF 
bonds. We have found that dX J — 0.125t the C-AF phase 
has higher energy than the CE one for small clusters 
of eight sites, and the same holds true for larger 4x4 
clusters. However, it has been pointed out that finite size 
effects are important;^ so one would like to investigate 
still larger clusters, at least at T = 0. 

As a first step, we exploit the quasi-lD nature of the 
two phases and investigate 8x2 clusters instead of 4 x 4 
ones (at present, we cannot treat more than 16 sites): (i) 
a ladder for the C-AF phase, and (ii) the cluster shown 
in Fig. [7{c) for the CE phase. While the energy of the 
CE phase hardly changes with cluster topology, the one 
of the C-AF phase is considerably lowered, albeit still 
higher than that of CE phase, see Tab. ^ To make the 
ground state calculations more conclusive and to elimi- 
nate systematic errors, it is therefore indeed necessary to 
investigate finite size effects. 

Due to the magnetic order in both C-AF and CE 
phases, hopping occurs only along a ID path, so chains 
can be investigated instead of 2D lattices. In this way 
larger systems could be reached. Taking the kinetic en- 
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ergy for a FM chain, 

^c/CE _ .i^y-ro-t ~ , ~t ~ 

i 

we considered two different geometries: (i) the C-AF 
phase with a chain along b axis, i.e., taking + sign for 
the interorbital hopping in Eq. 1321 and (ii) a zig-zag 
path chosen instead for the CE phase, which leads to a 
sequence (&, 6, a,a,b,b, . . .) of bond directions and there- 
fore to the phase sequence (+, +,—,—,+,+,...) in the 
hopping. 

Because the FM order within the chains (in either CE 
or C-AF phase) is perfect at zero temperature (u^ =1), 
the superexchange along them contains only one term. 



ttFM 



^E 



c^c 



2T^T' 
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'■i+l 



(33) 



where C is the directional orbital along the bond direc- 
tion, i.e., along b in the C-AF phase and alternating be- 
tween a and b in the CE phase. However, superexchange 
via AF bonds also contributes, and these interchain cou- 
pling terms could be crucial for stabilizing one or the 
other phase^S We therefore embedded the chains and in- 
cluded into the effective ID Hamiltonian additional AF 
superexchange terms. 



H^ = 



- i^niUi+i 



i 



the Hamiltonian comprising the three terms given in Eqs. 
(|S^ . (|S5|l . and (|S3Il on chains of different length. 

The results for chains of various lengths were next 
extrapolated to the 1/L -^ limit (see Fig. I10|) . An 
excellent fit was obtained using a quadratic dependence 
of the ground state energy on the inverse chain length, 
Eq = Eq^oo + k ■ yj. In this way we deduced the extrap- 
olated energy values, -Eq.oo- As only directional orbitals 
oriented along the chain contribute to the kinetic energy 
in the C-AF phase, this energy is not influenced by U. 
Consequently, one finds that large on-site repulsion U, 
i.e., small J, favors the C-AF phase (see also Refs. fis 49). 
Furthermore, it appears that the energies of both phases 
are so close to each other for J ^ t/8 that one cannot 
distinguish between these phases and decide on the na- 
ture of magnetic correlations in the ground state. This 
result is not strongly affected by a uniform crystal field 
Ez either — the energies of the two magnetic phases are 
again very similar for J = t/8. 

As the commonly used picture of the CE phase as- 
sumes charge order, one expects that it to be stabilized 
by nearest neighbor Coulomb repulsion V. The extrapo- 
lation to 1/L ^ gives indeed lower energies of the CE 
phase for F > 0, but the effect of V remains surprisingly 
small. The reason is that the second of the AF terms 
in Eq. (|34ll already induces some charge order in the 
C-AF phase as well, and therefore its energy does not 
suffer much from Coulomb repulsion V . If, on the other 
hand, V becomes too large, it hinders electron motion 
along the FM chains in both the C-AF and CE phases 
and thus affects both of them, even favoring the C-type 
AF phase for very large V > 1.5t. 
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J > n.-^n. 



iC'S+lC ' 



(34) 



where C is the directional orbital perpendicular to the 
chosen path, i.e., couples nearest neighbor sites on two 
adjacent chains. The form of Eq. (|^ is motivated by 
the fact that a bridge site on one chain (in CE phase) lies 
always next to two corner sites of the neighboring chains 
in ab plane. By symmetry, these should be equivalent to 
the corner sites on the considered chain, which are again 
the nearest neighbor sites to a given bridge site. The 
energy of the C-AF phase is evaluated with a similar 
term. 

In order to check the validity of this one-dimensional 
approach, we compare the energies obtained on small 2D 
clusters (third to fifth column in Tab.Hl) to those obtained 
on short chains of four and eight sites (last two columns 
in Tab.Hl). For the CE phase, L = 4 corresponds to the 
chain length encountered in a \/8 x \/8 cluster, and the 
energies differ indeed only by ~ 0.0015t. For the C-AF 
phase, L = 4 corresponds to either \/8 x -\/8 or 4 x 4 clus- 
ter, and one also finds remarkably good agreement. The 
difference is larger for L = 8 and the ladderlike clusters, 
but it seems still reasonable to investigate systematically 



C. Orbital and charge order at half doping 

When the magnetic order of CE type occurs, the sites 
along each FM zig-zag chain are nonequivalent, and one 
expects that holes are predominantly found at the corner 
gj|.ggj22j23j24 However, in the present finite cluster calcu- 
lations different CE patterns mix with each other and 
the holes cannot be detected using just density opera- 



TABLE I: Ground state energies as obtained for C-AF and 
CE phases with different clusters of Fig. Qand for ID chains 
simulating these phases, for three representative values of Ez . 
In case of 2 x 8 clusters, a ladder was used for the C-AF 
phase and the cluster shown in Fig. |7[c) for the CE phase. 
Parameters: J = 0.125t, V = 0. 
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-0.6466 
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FIG. 10: (Color online) Finite size extrapolation of the ground 
state energy Eq for the C-AF and CE phases on chains of 
different length L < 16, as obtained at T = for different 
values of J: (a) J = 0.25t, (b) J = 0.125t, and (c) J = 0.05i. 
The data for the C-AF phase alternate between the chains of 
length L — An and L — in + 2, with n integer. Parameters: 
J' ^E^ = V = 0. 



tors. This information can be extracted only from inter- 
site correlations, such as for instance the spin-hole-spin 
correlation TZ which measures the spin-spin correlations 
across a central site occupied by the hole, see Eq. (|29|l . 
For the clusters shown in Fig.[7|we consider the FM zig- 
zag chains along (11) direction. On the one hand, the 
sites to the 'right' and 'left' of a corner site along a axis, 
i.e., (10) direction, have opposite spin, as well as those 
'above' and 'below' it along b axis, i.e. (01) direction. 
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FIG. 11: (Color online) Orbital structure for the CE phase at 
T = 0, with circles (crosses) for 3x^ — r^ (3y^ — r'^) orbitals, as 
obtained for: (a) E^ = 0, (b) E^ = 0.2i, (c) E^ = 0.5t, and (d) 
Ez = t. The size of circles and crosses is proportional to the 
electron density in a given orbital. Parameters: J = 0.125t, 
y = 0. 



For the bridge sites, on the other hand, the neighboring 
spins along either direction have the same sign. In the 
CE phase, negative values of TZ indicate therefore that 
holes occupy corner rather than bridge sites. ^" 

We have calculated the correlation function 7?. for the 
ground state, i.e., first the spins were set into the zig- 
zag CE pattern and next the resulting orbital Hamilto- 
nian was solved with Lanczos diagonalization. In the 
absence of the nearest neighbor Coulomb repulsion (at 
V = 0), we then found TZ = —0.087, from which we de- 
duced the electron density at corner sites ric — 0.413 
vs. rifc ~ 0.587 on the bridge sites. (The MC data 
at low temperature f3t — 100 with J' = 0.05t and 
Ez = V = give a somewhat weaker spin/charge order 
with TZ = -0.0757 ± 0.001.) The electrons at the bridge 
sites are almost exclusively in the directional 3x^ — r^ 
{3y^ — r^) orbitals along a (6) axis, i.e., along the di- 
rection of the zig-zag chain. In contrast, the electrons at 
corner sites are more evenly distributed over both orbitals 
in the CE ground state (at i^z = 0), with n^ = 0.2307 in 
|a:) vs. n^ = 0.1823 in \z) orbital, respectively. 

These findings, i.e., the relatively small difference in 
the density at bridge and corner sites and the occupation 
of the directional orbitals, are in contrast to some experi- 
mental resultsi^iiS However, the 00 similar to our find- 
ings was also reported jiSiS*^ and one expects that rather 
extreme charge modulation, with holes at corner sites 
and eg electrons in bridge positions, should be excluded 
as then the FM double exchange which stabilizes the CE 
phase would be lost. Therefore, the charge order with 
alternating Mn'^^ and Mn^+ ions has been challenged, 
and intermediate valence picture with only small density 
variations has been suggested both in experimental^ and 
in theoretical studies i^^i^SiSiS 

In this Section we would like to address as well the 
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FIG. 12: (Color online) Possible orbital and spin structure 
for the CE phase, as suggested: (a) in Fig. 1 of Ref. 1 53 . and 
(b) from our calculations. The orientation of the rectangles 
indicates the type of occupied orbital - 
vertical y^ — z^ . 



recent controversy concerning the type of the 00 in the 
CE phase. We analyzed the orbital occupation in the 



Another difference between the analysis performed in 
Refs. |52 and our results is that charge order is not per- 
fect in our case, as it was assumed in their analysis. This 
does not influence the 00, however — iav V = t the 
charge order is enhanced (i.e. the electron density at 
bridge sites increases) without greatly affecting the type 
of the 00 (not shown). Here we would like to emphasize 
that our findings concerning the nature of spin and 00, 
as well as rather weak charge order, agree with the recent 
Hartree-Fock calculations on the multiband d-p model;^ 
indicating that the local Coulomb interactions and su- 
perexchange suffice already to stabilize the CE phase. 
There is no doubt, however, that the oxygens distortions 
also contribute to the stability of this phase^^*^ and 
would expand the regions of the CE phase in the phase 
diagrams shown in the next Section. At the same time 
realistic oxygens distortions would also modify somewhat 
the type of occupied orbitals, but the essential features of 
the orbital order in the CE phase, see Fig. ll2r b). would 
remain the same. 



3a;2 - r^ and Sy^ - r^ orbitals for various values of the ^^^ ^^^^^ *^^^ z -x jy -z - and 3x - t I'^V 



crystal field splitting i^^ > (Fig. [TTl)- At i;^ = 
the electrons at the bridge sites are found in the orbitals 
parallel to the FM chain [Fig. Ilir a)]. i.e., in 3x^ — r^ 
or 3y^ — T^ ^ as discussed above. One finds that the cor- 
responding orthogonal orbital is practically empty, e.g. 
2/^ — z^ orbitals at the bridge sites with both neighboring 
FM bonds along x direction. However, large density is 
found within these orbitals at the remaining bridge po- 
sitions (in 2/^ — z^ orbitals on those sites where the FM 
chains have the bonds along y direction), which is due to 
the considerable overlap between the 3j/^ — r^ and y^ — z^ 



orbitals. Thus, the ix^ - r"^ /iy"^ 



order appears to 



type orbital order are qualitatively similarjS and that 
one can come from one to the other one by adding a con- 
stant cn^tal field, are consistent with the results reported 
in Ref.lSll where a shear- type distortion (alternating con- 
tractions along x/y-axes) has been found more plausible 
than a JT type 00 (alternating elongations along x/y- 
axes). In the shear-type order, the out-of-plane Mn-0 
bond is of similar length as the longer in-plane bond, 
while it is comparable to the shorter bond in the JT-like 
case. Variation of the out-of-plane bond, equivalent to 
i?2 > in our model, can therefore lead from one sce- 
nario to the other. In closing, we remark that the kinetic 



be qualitatively similar to the z — x jy 



order. In energy favors ix^ 



V3y^ 



-occupation as this max- 



the present case, however, we would rather identify it 



as 3x2 _ ^2/3^2 



order because the occupation on 



the bridge sites is almost exclusively in the directional 
orbitals. 

When the crystal field E^ increases, the orbital occu- 
pation changes (see Fig. Ill|l . Not surprisingly, density 
in the iz^ — r^ orbitals (not shown) has increased, which 
modifies the type of the 00 — one finds now that the di- 
rectional orbital is the one which is empty on some sites, 
which means that the electron on a bridge site has moved 



imizes the hopping. Onsite Coulomb repulsion inhibits 
the kinetic energy and when we weaken Coulomb repul- 
sion from [/ = 8t (J = 0.125t) to C/ = 4i (J = 0.25i) 
and thus enhance the kinetic energy, we indeed find that 



stronger E^ is needed to induce z — x jy 



-type 00. 



from the Zx^ 



(3y — r^) into the z — x {y — z ) 



orbital. This situation could therefore be considered as 
representing z"^ — x^ jy^ — z^-type order. The observed 
transition from the former (at E^ = 0) to the latter (at 
E^. = V) phase is driven by the crystal field and is grad- 
ual, see Fig. ^2 Therefore, it may be argued that the 
combination of spin and orbital structure in Fig. 1 of Ref. 
152 . given schematically in Fig. ll2r a). is incorrect: In this 
picture, the FM spin chains run perpendicular to the oc- 
cupied z^ — x^ (y^ — z^) orbitals at the bridge sites, i.e., in 
y-[x-) direction. Instead, our data indicate that the FM 
chains run as shown in Fig. I12r b). i.e., in x-direction, 
where the z^ — x^ orbital is occupied for large E^ > Q 
and along y-direction for the y'^ — z^ sites. 



This is in accordance with the usual experience based on 
comparing LDA with LDA-I-U calculationsj^ where in- 
clusion of onsite interaction has likewise been found to 
be crucial in stabilizing shear-type over JT type 00. 



D. Phase diagrams at half doping 

Our results on the stability of magnetic phase &\, x — 
0.5 were collected in the phase diagrams of Fig. E|for two 
representative values of the on-site Coulomb repulsion, 
U = 8t and At, leading to J = t/8 and t/4. For a more 
realistic J — t/8, the extrapolated energies for the C-AF 
and CE phases are indeed very close to each other, with 
V > slightly favoring the CE phase (see above) , and we 
anticipate that the observed differences might actually be 
smaller than the error induced by the use of ID chains 
instead of 2D clusters. We therefore could not determine 
a phase boundary between these two phases. We also 
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FIG. 13: Phase diagram at half doping {x = 0.5) for: (a) 
J — 0.125f, and (b) J = — .25t, as obtained from finite size 
considerations at orbital degeneracy {E^ = 0) . The black lines 
around the phase boundaries give their estimated numerical 
and finite-size errors. At J/t = 1/8 the energies for the C- 
AF and CE phases differ only very little and we could not 
determine a phase boundary between them. For J/t = 1/4 
the CE phase has lower energy than the C-AF phase. 



had to compare the energies of the two quasi-lD phases 
to the 2D-like FM and AF phases. However, we observed 
that the energy of the CE phase for L = 8 is already very 
similar to the extrapolated result for L —f oo for all values 
of V. Therefore the results obtained for the 8x2 cluster, 
see Fig. EIc), could be used for energy comparison with 
the FM and AF phases. 

In order to estimate the uncertainty of the ground- 
state energies in a reliable way, and thus of the de- 
termined phase boundaries, we further investigated the 
change of the energy of the FM phase with the cluster 
size for three clusters used in the present calculations: 
VS X v^, VTO X yiO, and 4x4. It was found that the 
energies did not depend much on cluster size, but these 
might, however, still be too small. A somewhat simpler 
situation occurs for the G-type AF phase which turns out 
to be perfectly charge ordered for T^ > — J, so its energy 
is independent of V . Moreover, as all electrons are per- 
fectly localized, one finds in the present approximation a 



classically ordered G-AF phase without finite-size effects. 

For the G-AF and CE phases, we investigated how 
much the energy of the CE phase for L = 8 differs from 
the extrapolated energies for either the G-AF or the CE 
phase at L —^ oo. We finally arrived at the estimated er- 
ror AJ' < ^^ « O.Olt and the resulting phase diagram 
for J/t = 1/8 given in Fig. HSJa). The CE or G-AF 
phases are stable in between the AF and FM phases, re- 
spectively. Increasing nearest-neighbor Coulomb repul- 
sion V > weakly suppresses the FM phase and fa- 
vors the AF phase. On the one hand, the FM phase 
is suppressed because it is stabilized by the kinetic en- 
ergy which decreases when the charge order is induced 
by finite V. The AF phase, on the other hand, is already 
charge ordered in the absence of explicit Coulomb repul- 
sion and is therefore not influenced by V. Altogether, the 
electrons redistribute at finite V also in the CE phase, so 
its range shrinks. 

For a smaller on-site Coulomb repulsion U — At, e.g. 
J/t — 1/4, the CE phase is favored over the G phase 
for all values of V, at least in the investigated parameter 
range — 0.5i < V < 2t. We used the same procedure as 
described above to determine the boundary toward the 
FM and AF phases, and arrived at the phase diagram 
depicted in Fig. IT^ b). One finds that the G-AF phase 
extends to a broader range of parameters, and the region 
of CE phase is reduced with increasing V. This contra- 
dicts the common belief that the CE phase is stabilized 
by the nearest neighbor Coulomb interaction. 



E. Magnetic and orbital correlations for increasing 
doping 

We complete this Section by a qualitative discussion 
of the changes in orbital occupation and magnetic cor- 
relations in a monolayer under increasing doping. For a 
4x4 cluster doped with one hole {x = 1/16 = 0.0625) we 
observed rather weak AF order in MC runs for J' = 0.02t 
and more pronounced AF order for J' = 0.05t, with 
Ez — V = 0. Unfortunately, we were not able to per- 
form MC simulations for more than one hole on 4 x 4 
clusters (x > 1/16 — 0.065) due to the increasing size 
of the Hilbert space. For one hole doped to a smaller 
VS X VS cluster (doping x = 1/8 = 0.125), we found 
weakly AF correlations for J' — 0.05f, while AF order 
had vanished for J' = 0.02t, which is a clear sign of 
the increasing importance of FM double exchange with 
increasing doping. The fast disappearance of AF order 
{x > 0.125 for J' — O.Q2t) agrees with experimental ob- 
servations, where AF order disappears at x > 0.115, and 
is replaced by short-range spin-glass type of order »i£ 

Magnetic correlations at higher doping are theoret- 
ically challenging, and various anisotropic magnetic 
phases {A-AF and G-AF) were obtained in the model 
discarding strong Coulomb repulsion»^ In the range of 
large doping x > 0.5, higher doping can be reached 
for Ndi_2;Sri_|_a:Mn04, while Lai_a;Sri+j;Mn04 can be 
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FIG. 14: (Color online) MC results for the spin-spin corre- 
lations S{r) Ij26|l for a monolayer, as obtained in \/8 x \/8 
cluster with two electrons (at high x = 0.75 doping) for two 
parameter sets: J' = 0.05i, V ^0, E^ =0 [x for (10) and + 
for (11) direction]; J' = 0.025i, V = t,E^= 0.25t [A for (10) 
and V for (11) direction]. Error bars are smaller than symbol 
sizes. Parameter: J = 0.125t, f3t = 100. 



doped only up to a; « 0.7. In Ndi_a;Sri+2;Mn04 samples 
the C-AF phase was observed^? for 0.75 < a; < 0.9; addi- 
tionally, a structural phase transition suggesting predom- 
inant occupation of directional orbitals along one axis 
was reported. We found it very encouraging that the 
same trends have been observed in the MC simulations 
for two electrons on a \/8 x -\/8 cluster and for four elec- 
trons on a 4 X 4 cluster (x = 3/4), as well as for three 
electrons on 4 x 4 cluster (x — 0.8125), for large enough 
J' = 0.025i and 0.05i. 

The spin correlations obtained from MC simulations 
of a \/8 X \/8 cluster with two electrons are shown in 
Fig. El for two parameter sets: J' = 0.05i, V = Ez = 0, 
and J' = 0.025i, V = t, E^ = 0.25t. For a larger value 
of J' — 0.05t we see the telltale signal of the C-type 
phase, the strongly negative signal at (11). The nearest- 
neighbor spin correlation in (10) direction nearly vanishes 
because the FM and AF signals from the two directions 
almost cancel each other. (Parallel tempering was em- 
ployed in this case to ensure good autocorrelation times.) 
For smaller J' = 0.025i, and V = t, E^, = 0.25t, the C- 
AF phase is not as marked, and the FM correlations are 
stronger than the AF ones. In agreement with expec- 
tations based on the ID model,-- the electrons occupy 
directional orbitals along the FM direction in the C-AF 
phase, because this maximizes the gain of the kinetic en- 
ergy. 

Figure 1151 shows that orbital occupation depends 
strongly on doping for realistic parameters, J' = 0.025i 
and V ^ t. This is reflected by the percentage of elec- 
trons occupying out-of-plane \z) orbitals which is further- 
more very sensitive to the actual value of E^, see Fig. El 
It is quite remarkable that for the degenerate Cg orbitals, 
i.e., without crystal field {Ez = 0), practically only in- 
plane \x) orbitals are occupied at a; = [squares in Fig. 
I15r a'l]. This state is induced by finite AF superexchange 
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FIG. 15: (Color online) Orbital electron densities for increas- 
ing doping X = 1 — n as obtained in MC method for a y/S x y/S 
cluster: (a) percentage of density in \z) orbitals riz/n for a few 
selected values of crystal field parameter Ez and for variable 
Ez (see below); (b) absolute densities riz and n^ in the two or- 
bitals from MC (symbols) compared with the densities found 
with the KR slave boson mean-field approximationSi (lines) 
with variable Ez = {\ —x)t. The charge distribution found at 
doping X = 0.5 with the KR method (n^ ~ 0.14, n^ — 0.36) 
is indicated by diamonds. Remaining parameters in MC cal- 
culations: J = 0.125f, J' = 0.025t, V ^t, pt^ 100. 



J' — 0.025i which (due to large Hund's exchange) selects 
the AF interactions between Cg, and these interactions 
are maximal in ab planes when \x) orbitals are occupied. 
In this way the core spin superexchange influences also 
the Cg orbital occupation (at J' = one finds an al- 
most isotropic electron distribution with Uz/n ~ 0.5 in 
the FM phase). Upon doping, electron population in 
\z) orbitals gradually increases, reaches a maximum at 
X — 5/8 where FM correlations are found, and then de- 
creases again. This behavior follows from the kinetic en- 
ergy which contributes in the entire regime of < a; < 1 , 
and is gained when the interorbital processes which ex- 
cite electrons to \z) orbitals, ^ \/3( '^ 
allowed. 
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For a large positive crystal field favoring \z) orbitals 
Ez — 0.5t [triangles in Fig. Iisr a)]. the density distri- 
bution is reversed at a; = — almost all electrons are 
found within \z) orbitals. They redistribute, however, 
gradually with increasing doping, because the kinetic en- 
ergy competes with the crystal field and favors in-plane 
\x) orbitals. (In this case one finds a FM phase in a range 
of doping from x = 1/4 to x = 5/8, separated by the CE 
phase found at a; = 0.5 for i?^ = 0.) In the intermedi- 
ate case of a smaller crystal field splitting of E^ = 0.2i 
(crosses), the normalized \z) density nz/n is first reduced 
but next rises again towards the FM phase which in this 
case suppresses the CE phase at a; = 0.5. 

As the last scenario, we investigated the effect of a 
doping-dependent crystal-field splitting. 



x]t. 



(35) 



The doping dependence of this type [circles in Fig. ll5r al] 
is qualitatively expected by considering the experimental 
datai^ This is probably the most realistic case of those 
considered here, as it shows both the correct orbital po- 
larization for X = and the CE phase for x = 0.5. Note 
that in all cases except E^ = 0, we find the experimen- 
tally observed*- increase of in-plane \x) electron density 
with increasing doping which shows that Ez > in the 
low doping regime. 

Taking the decreasing with x crystal field as in Eq. 
(|35|l . one finds a very fast decrease of n^ from n^ == 1 at 
X = to riz ~ 0.22 at a; = 0.25 [Fig. Elb)]. Thus the 
electrons move fast from \z) to [a;) orbitals in this doping 
regime, as dictated by the kinetic energy gain. This qual- 
itative trend is very well reproduced by the analytic ap- 
proach using the mean-field approximation in the slave- 
boson method, introduced by Kotliar and Ruckenstein.^^ 
We have adapted this method to the FM phase in a mono- 
layer, as explained in the Appendix. The performed com- 
parison with the results of exact diagonalization shows 
that this analytic approach provides a surprisingly reli- 
able way to estimate both the charge distribution in lay- 
ered systems and the magnetic interactions at increasing 
doping. In the present case the double exchange mech- 
anism promotes FM states in a large range of doping, 
while in a bilayer system changing electron density dis- 
tribution provides a natural explanation for an observed 
transition from the FM to A-AF structure;-^ 

The changes in magnetic correlations at increasing 
doping follow from the competing superexchange and 
double exchange interactions. The double exchange is 
directly proportional to the kinetic energyj24 which van- 
ishes at a; = and is gradually gained by doping when 
the available space for hopping processes increases up to 
X ~ 0.5, and then is gradually lost beyond half doping 
(see Fig. I16|) . Therefore, the kinetic energy has an ap- 
proximately parabolic shape, as obtained for a v8 x v8 
cluster with a few representative parameter sets. Since 
different orbital occupation can favor either FM or AF 
spin configuration for any bond, the superexchange en- 



UJ 




FIG. 16: (Color online) Kinetic (approximately parabolic 
shape) and superexchange (weakly increasing) energies for in- 
creasing doping X = 1 — n, as obtained for %/8 x \/8 clusters 
with various parameter sets: O s^nd O — for J' = 0, and 

£z = V = 0; X and H for J' = 0.05i, and E^ = V = 0; A 

and V — for J' = 0.025i, V ^ t, and E^ given by Eq. (ESJ. 
Parameters: J = 0.125t, pt = 100. 



ergy is finite even in the cases where the average nearest- 
neighbor spin correlation function vanishes, like in the 
CE and C-AF phases. 

In the FM case at J' = 0, the kinetic energy has its 
minimum at a; = 0.5 because of the optimal carrier den- 
sity fFig. 1161) . For J' > 0, however, the minimum of the 
kinetic energy is moved to larger doping x — 5/8 — 0.625, 
because this electronic filling allows to realize the FM 
phase, while the magnetic correlations favor instead the 
CE phase at a; = 0.5 at the expense of the kinetic energy. 
This demonstrates that these two energies are to some ex- 
tent complementary and their competition controls the 
magnetic order. With finite nearest neighbor Coulomb 
repulsion V = t, the kinetic energy of the CE phase is 
partly lost due to charge order. Concerning the total 
energy, which contains additional contributions from V 
and Ez and is not shown in Fig. 1151 we want to empha- 
size that it is convex over the entire doping range and 
for all parameter sets. Although this result seems to ex- 
clude phase separation, ^^ the clusters used in the present 
study are definitely too small to address this issue in a 
conclusive way. 



IV. MAGNETIC PHASES IN BILAYER 

MANGANITES 

A. Phase diagrams for undoped system 

Bilayer manganites like La2-2a;Sri_|_2a;Mn207 represent 
an intermediate situation between the 2D monolayer sys- 
tems and 3D perovskite manganites, so it is interesting to 
ask to what extent the qualitative trends reported above 
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FIG. 17: Phase diagram of the undoped bilayer system as 
obtained in {E^, J') plane with a -/S x \/8x 2 cluster at T = 0. 
The phases found between FM and G-AF order for two ab 
planes with interlayer coupling along c axis are: C-AF and 
A-AF phases as in Fig. Q] C'-AF: FM along c axis and AF 
within the ab planes; A'-AF: FM along a and c axes and AF 
along 6 axis. The 00 which accompanies the FM and C-AF 
phase at _Ez = is also indicated, and shown in Figs. Il^f a)- 
[THIc). Units J'/t and E^/t correspond to J = 0.125t. 



for the monolayers are modified by the interlayer cou- 
pling. We shall provide some limited answers to this 
question, as unfortunately the calculations could only be 
performed for certain selected fillings of the smallest bi- 
layer VS X -s/S X 2 cluster used in the numerical studies. 

The study of magnetic correlations in an undoped 
(x — 0) system performed on VS x a/8 x 2 clusters led to 
the phase diagram in the {J',Ez) plane (see Fig. I17|l . It 
was obtained by comparing the energies of various pos- 
sible magnetic states at T = 0. As in a monolayer, large 
positive J' leads to the G-AF phase (with AF correlations 
in all three spatial directions), and strongly negative J' 
induces FM phasen^ Apart from these two phases, one 
finds here several different types of intermediate magnetic 
order in the crossover regime from FM to G-AF phase, 
with FM bonds along either two or only one cubic direc- 
tion, see Fig. ^ Two phases of the same type, either A- 
AF and ^'-AF, or G-AF and G'-AF, are distinguished by 
the directions of FM bonds. When the FM correlations 
occur within the ab layers, as found experimentally ^AS, 
they are called ^-AF and G-AF phases, respectively. 
These phases appear for the expectedi^ positive {E^ > 0) 
crystal field which favors \z) orbital occupancy, while for 
Ez < interlayer FM correlations occur in A'-AF and 
G'-AF phases. 

For intermediate J', one obtains the ^-AF phase, rem- 
iniscent of the ground state of undoped 3D perovskite 
manganite LaMnOa. For Ez > 0.2t this plane is partic- 
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FIG. 18: (Color online) Orbital correlations Te{r) ^ in (01) 
direction within one ab plane and between the two different 
planes of the bilayer, as obtained for the undoped \/8 x \/8 x 2 
clusters at S^ = and T = in: (a) G-AF phase, (b) FM 
phase, and (c) A-AF phase. In each case Tf operators are 
defined by different bases of orthogonal orbitals: ^ = — 
{\x), \z)}, and 8 = n/2 — {|a;) -|- \z), \x) — \z)}. Parameters: 
(a) J' = 0.05t; (b) and (c) J' = -0.02t. 



ularly robust which suggests that positive Ez simulates 
here the effect of missing planes along c axis on the elec- 
tron distribution. Large \z) amplitude is needed as (i) 
this type of order is supported by the AO order in the ab 
planes, and (ii) \z) orbitals are responsible for the AF in- 
terlayer coupling. This coexisting AO/FM phase within 
ab planes is robust and can be suppressed only by AF 
core spin superexchange J' > of a similar value as in 
the monolayer (see Fig. |2J|. 

MC simulations performed at f3t — 100 for various pa- 
rameter sets support the phases found in the ground- 
state calculations (J = t/8 in all cases): J' = 0.05t, 
Ez ^ (G-AF); J' = -0.02t, Ez = (FM); J' = 
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and J' = -0.02i, E, = 0.5i (A-AF); J' = -0.005t, 
Ez = — 0.5t (C"-AF). For some parameter sets close to 
the ground-state phase boundaries, we observe compet- 
ing phases, i.e., configurations showing several different 
phases occur in the MC runs: J' — 0, E^ — (FM and 
A-AF), J' = O.Oli, S^ = {A-KF, some configurations 
with A'-AF and C-AF). 

The crossover region from the FM to the G-AF phase 
with increasing J' is characterized by the competition 
between nearly degenerate magnetic phases. As for 
the monolayer phase diagram, the ground-state calcu- 
lations also yield small regions with the C-AF or C'-AF 
phases. The C'-AF was indeed found in MC runs for 
J' = -0.005i, Ez = -0.5i, I3t = 100, J = l/8i, but the 
C-AF phase (having a very narrow stability region in the 
ground-state phase diagram) only occasionally surfaced 
in the MC runs for E^ = and J' = O.Olt (competing 
with A-AF) and J' = 0.02t (competing with G-AF). It 
also occasionally occurs at J' = 0.015t, but the mag- 
netic structure there remains unclear, which could mean 
that temperature was still too high and/or the cluster 
too small. 

For negative Ez in-plane \x) orbitals are favored, which 
in turn enhances the region of stability of the G-AF 
phase. This can be understood by looking at the orbital 
correlations in the G-AF phase at i^^ = depicted in 
Fig. Iisr al : Even at orbital degeneracy (without a crystal 
field) , one finds rather pronounced polarization for ^ = 
— 80% of the electrons occupy \x) orbitals as then the 
superexchange energy is gained. The G-AF phase can 
therefore take advantage of a negative E^ , as seen in the 
phase diagram of Fig. 1171 

Figure ITHT b) shows the orbital correlations in the FM 
phase, where strong AO order within the layers coex- 
ists with FM order in ab planes. As this AO order can 
best develop when occupations of \x) and \z) orbitals are 
nearly equal, the FM phase is most pronounced around 
Ez = 0. At the same time, the interlayer coupling (ki- 
netic energy) is much weaker. This state competes with 
the G-AF state with larger occupancy of \x) orbitals. 

Finally, the A-AF phase is found mainly for Ez > 
which enhances electron density in \z) orbitals and sup- 
ports the AF interlayer coupling. The orbital correlations 
in this phase (at Ez = 0) are depicted in Fig. ITHT c) — 
one finds the AO order within the FM planes, i.e., or- 
bital correlations for 6 = 7r/4 are strongly alternating 
within the ab planes. Between the planes, i.e., along the 
AF bonds in c-direction, the orbital correlation is weakly 
positive, indicating weak FO order. 

Although the phase diagram of Fig. 1171 was obtained 
with small V8 x -\/8 x 2 clusters, we argue that it is rep- 
resentative of the bilayer system in the thermodynamic 
limit, as it is determined by short-range spin and orbital 
correlations that follow from local superexchange interac- 
tions (charge excitations) on the nearest-neighbor bonds. 
To support this point of view we present also the phase 
diagram found with a smaller 2x2x2 cluster in Fig. 
[T^ One finds that indeed the same phases occur as in 
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FIG. 19: Phase diagram of the undoped 2x2x2 cluster 
as obtained in {E^, J') plane at T = 0. The phases and the 
parameters are the same as in the bilayer system (Fig. I17L 
with the directions of FM bonds in yl-AF, A'-AF, and C'-AF 
phases distinguished by the symmetry-breaking crystal field 
oc Ez, see Eq. Ijlljl . C-AF phase almost vanishes in this case. 



Fig. 1171 and their stability regimes are remarkably close 
to those of the larger VS x \/8 x 2 cluster. 



B. Competition betw^een different phases in 
half-doped bilayer clusters 

Next, we investigate the magnetic and orbital order at 
half doping, assuming orbital degeneracy {Ez = 0). Fig- 
ure I^OIa) shows the energy of various phases of the half- 
doped bilayer depending on the value of the f 2g superex- 
change J' in the absence of intersite Coulomb repulsion 
{V — 0). Not only for negative but also for small positive 
J', the system is FM, while for larger positive J' > 0.02i 
one finds G-AF phase. In between these phases the CE 
phase, with alternating FM zig-zag chains in ab planes 
and AF coupling between them, has a lower energy. How- 
ever, finite-size effects are again important, as discussed 
in Sec. 



nil Bl and thus one expects that instead the C-AF 
phase is the actual ground state in the thermodynamic 
limit for the present realistic parameters. (Similar to the 
ID chains used for finite-size considerations instead of 
2D clusters, one can use 2D ladderlike clusters instead of 
3D clusters, but as the attainable lengths are here much 
shorter definite results are difficult to obtain. However, 
energy comparison of the C-AF and CE phases on 4 x 2 
and 8x2, including 6x2 for C-AF phase, suggests that 
the C-AF phase has lower energy.) 

We have already shown in the case of a monolayer that 
the nearest neighbor Coulomb repulsion V decreases the 
range of stability of the CE phase. Also for the bilayer 
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FIG. 20: (Color online) Ground-state energy E of various 
magnetic phases for increasing t2g superexchange J' in the 
half-doped bilayer \/8 x \/8 x 2 cluster, as obtained at T = 
for: (a) V — 0, and (b) V — t. Different phases with 
coexisting AF and FM bonds are defined as follows: C-AF 
— FM in a direction; A-AF — AF in c directions, FM in 
ab planes; CEi — FM zig-zag chains in ab planes with FM 
interlayer coupling; CE2 — FM zig-zag chains in ab planes 
with AF interlayer coupling. Parameters: J — 0.125i and 
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FIG. 21: (Color online) Spin correlations 1261 as obtained 
from MC simulations with a y/S x \/8 x 2 cluster in the highly 
doped regime for increasing coordinate r in ab plane: (a) G- 
AF phase for one electron (x = 15/16), and (b) C-AF phase 
for three electrons. Different symbols and lines indicate dif- 
ferent types of neighbors: x and dashed (blue) — within 
the ab planes along either (01) or (10) direction; o and solid 
(red) — within the ab planes along the (11) direction; D and 
dashed (green) — between the two planes along a or 6 di- 
rection; A and solid (black) — between the two planes and 
along the (11) direction in ab plane. Parameters: J — 0.125f, 
J' = 0.05f, V = E:,=0, and /3i = 100. 



system the CE phase is suppressed hy V ^ t, with weak 
preference of the C-AF phase [Fig. I2()r b)]. Namely, one 
finds that the energy of the C-AF phase is lower than that 
of the CE phase for ^ = i and J' > 0.02t. The reason is 
that the interlayer charge stacking (along c direction) re- 
quired by the CE phase costs extra energy now, whereas 
the C-AF phase permits alternating charge order in all 
three directions and has thus a lower energy-increment 
due to V. Unlike for the 2D clusters of Sec. IIII Dl the 
situation here is also similar for J = 0.25i (not shown). 
We have verified that both phases here have very similar 
energies on an 8 x 2 ladder, which seems to indicate that 
the C-type phase wins in the thermodynamic limit for 
J = i/4 as well. For ^ = i and J' < 0.02t the A-AF 
phase is more stable than the CE phase, which suggests 
that this parameter range is relevant for the experimen- 
tally measured LaSr2Mn2 07 sample j^^i^^ 

However, isotropic magnetic phases, FM for J' < 



0.02i and G-AF for J' > 0.02t, are more stable than 
anisotropic phases in the entire range of J', if V — t 
[Fig. I20r b)]. As for monolayers, finite ^ > favors the 
C-AF phase, because its kinetic energy vanishes already 
for y = and thus the energy of this phase is not af- 
fected by V , while charge order induced by V in other 
phases hinders electron motion, and thus their total en- 
ergies increase. While the spin-orbital model leads to the 
CE phase in a monolayer at J > t/8 or for y > 0, lattice 
degrees of freedom are apparently needed to stabilize it in 
3D perovskites. Furthermore, it should be noted that the 
bilayer compound LaSr2Mn207 shows A-AF order jiiii^ 
in contrast to the monolayers and the 3D compounds. 
The origin of this behavior remains unclear at present, 
especially as the orbital correlations are reported to be 
similar in all cases liSsi^ 
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C. Bilayer clusters at large doping 

The reference system for large hole doping regime is 
X = 1 case, when itinerant eg electrons are absent and the 
magnetic order depends exclusively on core spin superex- 
change J' — one finds then the G-AF phase for J' > 
(while an unphysical J' < induces the FM phase). We 
have been able to investigate the highly doped regime 
because the small number of electrons leads here to a 
small Hilbert space which allows to perform MC simula- 
tions down to a; > 3/4 on a -\/8 x VS x 2 cluster (filled 
by up to 4 electrons). MC simulations show that G-AF 
order for J' = 0.05i is strong enough to persist upon 
inclusion of one eg electron, i.e., at doping x — 15/16, 
as demonstrated by spin correlations presented in Fig. 
I2ir a'). They describe G-type alternating spin order in all 
three directions. 

Somewhat lower doping x — 13/16 (three electrons in 
the present cluster) gives a G-AF phase with FM chains 
lying in the ab planes, and predominant occupation of 
directional orbitals along the FM direction — this state 
is stabilized by the double exchange mechanism. Indeed, 
the strongly negative spin correlation at (1,1) point, see 
Fig. I21f b). is a signature of the G-AF phase. For distance 
r = 1 along (0, 1) direction, the signal is approximately 
zero, because the FM and the AF correlations in the two 
directions nearly cancel each other. Finally, the inter- 
mediate case of a; = 0.875 (with two electrons) does not 
allow to establish a clear picture of magnetic correlations 
(not shown), which may be due to strong competition 
between the two AF phases. In fact, experiments show 
a transition from G-AF to G-AF at doping x « 0.9 (see 
Ref. 12') which agrees with these results. 



V. DISCUSSION AND CONCLUSIONS 

The present study clarifies that orbital degrees of 
freedom are of crucial importance for the understand- 
ing of magnetic correlations in layered manganites. We 
treated a realistic model including intra- and interorbital 
Coulomb interactions and investigated charge, intcrsitc 
spin and intersite orbital correlations in monolayer and in 
bilayer manganites. The obtained results revealed a close 
relationship between orbital and magnetic order which 
follows the Goodenough-Kanamori rules at cc = 04^ The 
magnetic phases found in different doping regimes, where 
double exchange also contributes, are in accordance with 
experiments over the whole doping range < x < 1, 
particularly for the monolayer systems. 

For the undoped monolayers the model predicts either 
FM or AF order, but not the E-AF phase, reported previ- 
ously in an approach similar to ours but ignoring on-site 
Coulomb repulsion between the eg electrons i22i2^ In their 
study it was stabilized by the kinetic energy and arose 
mainly for nearly vanishing electron-phonon coupling,-"^ 
i.e., in the situation when lattice degrees of freedom could 
be neglected. Otherwise, the model of Ref. |23 is similar 



to our Hamiltonian apart from missing Coulomb repul- 
sion. Experimentally, however, the E-AF phase is only 
observed^ for the very strongly JT distorted HoMnOs, 
and never in less distorted compounds. In the FM phase, 
the 00 induced by the orbital superexchange, i.e., by 
local Coulomb repulsion, is found even in the absence 
of electron-phonon coupling, in contrast to the model 
without Coulomb repulsion. ^■^ This shows that the cor- 
rect treatment of electron correlation effects due to large 
Coulomb repulsion, which suppresses the kinetic energy 
in undoped compounds (at x = 0), is crucial for the qual- 
itatively correct description in this doping regime. 

The experimental situation in doped 

Lai_2,Sri+2,Mn04 could be modeled with varying 
crystal field favoring out-of-plane \z} orbitals in the 
undoped system, and gradually decreasing with x to 
accelerate the electron transfer from \z) to |a;) orbitals. 
Indeed, for positive E^ '^ 0.5i the undoped monolayers 
contain then almost only \z) electrons, while \x) occu- 
pation grows rapidly with doping when Ez decreases in 
the present model, see Fig. ^| in Sec. IIII El Indeed, 
such a doping dependence of E^ is suggested by recent 
experiments]^ 

Another success of the model is that one observes the 
CE phase at half doping with physically realistic param- 
eters for layered manganites, i.e., it is obtained for small 
t2g superexchange J' > 0.03i, as deduced^^ from the 
analysis of exchange constants in LaMnOs. We also in- 
vestigated the impact of nearest neighbor Coulomb repul- 
sion and found it to slightly stabilize the CE phase with 
respect to the G-AF phase in monolayers in the relevant 
regime of J', but to favor instead the G-AF phase in bi- 
layer clusters. (In both cases, the stability region of either 
the CE or the G-AF phase shrinks and that of the G-AF 
phase grows when nearest neighbor Coulomb repulsion 
is included.) In the CE phase, we found relatively simi- 
lar electron densities at corner (n^ ~ 0.413) and bridge 
{n^ ~ 0.587) positions in the zig-zag FM chains, which 
clearly contradicts the localized picture of this phase. 
Furthermore, we observed that electrons at the bridge 
sites are found in the directional 3x^ — r^/Sy^ — r^ or- 
bitals without crystal field and the planar z^ — x^ jy^ — z^ 
orbitals for Ez > 0- While this is in contrast to the in- 
terpretation of some experiments, ^^'^^ other groups re- 
ported similar charge distribution^ It is quite remark- 
able, however, that the CE phase could be here explained 
by a purely electronic mechanism — one may expect that 
the oxygen distortions due to the JT effect would further 
stabilize it. Also at large doping x ~ 0.75, the calcula- 
tions for monolayer clusters predict the G-AF phase with 
predominant occupation of directional orbitals, in agree- 
ment with experimental data for Ndi_a;Sri+2,Mn04iffi 

An interesting variation of spin and orbital correla- 
tions with doping was also found in the bilayer systems. 
They can be considered as intermediate between 2D and 
3D manganites, and we obtained the A-AF phase for the 
realistic parameters at x = 0, observed in undoped 3D 
LaMnOa perovskite compound. The absence of the CE 
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phase in the bilayer phase diagram of Ling et alm^ could 
not be explained, however. Perhaps the electron transfer 
from \z) to \x) orbitals at increasing doping is really fast, 
as suggested by the variation of intralayer and interlayer 
exchange constants, ^^■"'^^ and then the ^-AF phase is sta- 
bilized again, yet by different physical mechanisms. For 
very large doping x > 0.75, however, we obtained the 
C-AF and G-AF phases with a transition between them, 
as indeed observed in bilayer compounds .i^ 

Summarizing, the present study shows that the inter- 
nal frustration of magnetic interactions in doped man- 
ganites, with competing FM/AF terms in the spin su- 
perexchange which coexist with complementary terms in 
the orbital superexchange, has important consequences. 
Due to the intricate energy balance between different 
types of intersite correlations, the magnetic order may be 
completely switched over by small changes of microscopic 
parameters, when the orbital order which coexists with it 
switches at the same time. We find it quite encouraging 
that these generic features, as well as the experimentally 
observed trends in layered manganites, could be repro- 
duced within the present microscopic model. 
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where c]^ are the corresponding creation operators, and 
Xa — +27r/3, Xb = — 27r/3 are the phase factors for the 
bonds (ij) along a and b axis. Note that the crystal field 
term contains only off-diagonal terms for the complex 
orbital stares (|A.1|I . In order to implement rigorously 
the constraint at t/ = cx) we replace now the fermion 
operators as follows. 



4 = 44^.' (A-3) 



corresponding to a representation of the local states by 
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where |vac) is a true vacuum, following Ref. |2l|. In the 
slave-boson mean-field approximation we replace the bo- 
son operators by their averages, which leads to the (a pri- 
ori site-dependent) hopping rcnormalization factors qi±. 
For an isotropic charge distribution one finds then from 
a global constraint. 
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APPENDIX: SLAVE BOSON APPROACH FOR 
THE 2D MODEL OF SPINLESS FERMIONS 

It is notoriously difficult to implement electron correla- 
tion effects in nonmagnetic phases realized in multiband 
models. Therefore we consider here a simpler case of a 
FM monolayer in the limit of large [/ ^ oo to compare 
the resulting charge distribution with the exact diago- 
nalization of finite 2D clusters. The electronic structure 
for eg electrons is then described by the so-called orbital 
Hubbard model of Ref. ill 

It was shown recentlyiSi that cubic invariance is obeyed 
when the constraint of no double occupancy in the limit 
of large U is implemented by slave bosons for electronic 
states, using a basis of complex orbitals at each site i, 

|+). = 73(k). - ^\X)^), \-)^ = -^{\Z)^ + i\x),). 

(A.l) 
Then the hopping term Hf^ H14|) and the crystal-field 
term Hz Hll|) may be written as follows (the superex- 
change terms vanish in the limit of C/ — > oo): 
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are the same for all kinetic energy terms. In this way one 
arrives at the effective Hamiltonian with renormalized 
hopping terms. 
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(A.7) 



By diagonalizing it in reciprocal space and using the in- 
verse transformation to Eq. (jA.ip . 



= M 



= M 



(A., 



we determined the occupations of \x) and \z) orbitals 
shown in Fig. I15f b). A similar analysis used before 
for the bilayer system gave the density distribution and 
the effective exchange constants in good agreement with 
experiment li^ 
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